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Abstract 

An advantageous feature of piecewise constant policy timestepping for Hamilton-Jacobi- 
Bellman (HJB) equations is that different linear approximation schemes, and indeed different 
meshes, can be used for the resulting linear equations for different control parameters. Standard 
convergence analysis suggests that monotone (i.e., linear) interpolation must be used to transfer 
data between meshes. Using the equivalence to a switching system and an adaptation of the 
usual arguments based on consistency, stability and monotonicity, we show that if limited, 
potentially higher order interpolation is used for the mesh transfer, convergence is guaranteed. 

We provide numerical tests for the mean-variance optimal investment problem and the uncertain 
volatility option pricing model, and compare the results to published test cases. 

Key words: fully nonlinear PDEs, monotone approximation schemes, piecewise constant policy 
time stepping, viscosity solutions, uncertain volatility model, mean variance 
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1 Introduction 


This article is concerned with the numerical approximation of fully nonlinear second order partial 
differential equations of the form 


0 = F{x,V,DV,D^V) 


Vr - supqgQ LgV, X e M"* X (0, T], 
F(x) - g(x), X G X {0}, 


( 1 . 1 ) 


where x = (S, r) contains both ‘spatial’ coordinates S € and backwards time r. For fixed qm. a, 
control set Q, Lq is the linear differential operator 

LqV = tx{aqalD‘^V) +fi'^DV- VqV + fq, (1.2) 

where Uq G fj,q G M'’* and r^, G M are functions of the control as well as possibly of x. An 

initial (in backwards time) condition F(0, •) = ^(•) is also specified. 
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These equations arise naturally from stochastic optimization problems. By dynamic program¬ 


ming, the value function satisfies an HJB equation of the form (1.1). Since dynamic programming 


works backwards in time from a terminal time T to today t = 0, it is conventional to write PDE 


(1.1) in terms of backwards time t = T — t, with T being the terminal time, and t being forward 


time. 


Many examples of equations of the type (1.1) are found in financial mathematics, including 


the following: optimal investment problems 132]; transaction cost problems HZ); optimal trade 
execntion problems [T]; values of American options [25|; models for financial derivatives under 
uncertain volatilities mEoi; utility indifference pricing of financial derivatives m- More recently, 
enhanced oversight of the financial system has resulted in reporting reqnirements which include 
Credit Value Adjustment (CVA) and Funding Valne Adjustment (FVA), which lead to nonlinear 


control problems of the form (1.1) [I2ll3llll3|. 

If the solution has sufficient regularity, specifically for Cordes coefficients, it has recently been 
demonstrated that higher order discontinuous Galerkin solutions are possible |36|. Generally, how¬ 
ever, these problems have solutions only in the viscosity sense of m- 

A general framework for the convergence analysis of discretization schemes for strongly nonlinear 


degenerate elliptic equations of type (1.1) is introdnced in [7|, and has since been refined to give error 
bounds and convergence orders, see, e.g., laisiEi. The key requirements that ensure convergence 
are consistency, stability and monotonicity of the discretization. 

The standard approach to solve (1.1) by finite difference schemes is to “discretize, then opti¬ 


mize”, i.e., to discretize the derivatives in (1.2) and to solve the resulting finite-dimensional control 


problem. The nonlinear discretized equations are then often solved using variants of policy itera¬ 
tion m, also known as Howard’s algorithm and equivalent to Newton’s iteration under common 
conditions [9|. 

At each step of policy iteration, it is necessary to find the globally optimal policy (control) at 
each computational node. The PDF coefficients may be snfhciently complicated functions of the 
control variable q such that the global optimum cannot be found either analytically or by standard 
optimization algorithms. Then, often the only way to gnarantee convergence of the algorithm is to 
discretize the admissible control set and determine the optimal control at each node by exhanstive 
search, i.e., Q is approximated by finite subset Qh = {Qi,---qj} C Q. This step is the most 
computationally time intensive part of the entire algorithm. Gonvergence to the exact solution is 
obtained by refining Qh- 

Of course, in many practical problems, the admissible set is known to be of bang-bang type, 
i.e., the optimal controls are a finite subset of the admissible set. Then the trne admissible set is 
already a discrete set of the form Qh- 

In both cases, if we use backward Fuler timestepping, an approximation to at time 

is obtained from 


yn-i-l _ yr- 

At 


— max L„ 
qj&Qn 


= 0 , 


(1.3) 


where we have a spatial discretization L ^,, with h a mesh size and Ar the timestep. 


1.1 Objectives 


It is our experience that many industrial practitioners find it difficult and time consnming to 
implement a solution of equation (1.3). As pointed out in [3l|, many plausible discretization 
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schemes for HJB equations can generate incorrect solutions. Ensuring that the discrete equations 
are monotone, especially if accurate central differencing as much as possible schemes are used, is 
non-trivial |38j . Policy iteration is known to converge when the underlying discretization operator 
for a fixed control is monotone (i.e., an M-matrix) [9]. Seemingly innocent approximations may 
violate the M-matrix condition, and cause the policy iteration to fail to converge. 

A convergent iterative scheme for a finite element approximation with quasi-optimal convergence 
rate to the solution of a strictly elliptic switching system is proposed and analysed in [lOj . Here, we 
are concerned with parabolic equations and exploit the fact that approximations of the continuous¬ 
time control processes by those piecewise constant in time and attaining only a discrete set of 
values, lead to accurate approximations of the value function. 

A technique which seems to be not commonly used (at least in the finance community) is based 
on piecewise constant policy time stepping (PCPT) [28l|6]. In this method, given a discrete control 
set Qh = {(?i, • • • qj}, J independent PDEs are solved at each timestep. Each of the J PDEs has 
a constant control qj. At the end of the timestep, the maximum value at each computational node 
is determined, and this value is the initial value for all J PDEs at the next step. 

Convergence of an approximation in the timestep has been analyzed in m using purely prob¬ 
abilistic techniques, which shows that under mild regularity assumptions a convergence order of 
1/6 in the timestep can be proven. In this and other works |26[ I28j . applications to fully discrete 
schemes are given and their convergence is deduced. These estimates seem somewhat pessimistic, 
in that we typically observe (experimentally) first order convergence. 

Note that this technique has the following advantages: 

• No policy iteration is required. 

• Each of the J PDEs has a constant policy, and hence it is straightforward to guarantee a 
monotone, unconditionally stable discretization. 

• Since the PCPT algorithm reduces the solution of a nonlinear HJB equation to the solution 
of a sequence of linear PDEs (at each timestep), followed by a simple max or min operation, 
it is straightforward to extend existing (linear) PDE software to handle the HJB case. 

• Each of the J PDEs can be advanced in time independently. Hence this algorithm is an ideal 
candidate for efficient parallel implementation. 

• In the case where we seek the solution of a Hamilton-Jacobi-Bellman-Isaacs (HJBI) PDE of 
the form 

Vr — inf sxvpLqpV = 0, (1-4) 

the discretize and optimize approach may fail due to the fact that policy iteration may not 
converge in this case IS?!. However, the PCPT technique can be easily applied to these 
problems. 

In view of the advantages of piecewise policy time stepping, it is natural to consider some 
generalizations of the basic algorithm. Since each of the J independent PDE solves has a different 
control parameter, it is clearly advantageous to use a different mesh for each PDE solution. This 
may involve an interpolation operation between the meshes. If we restrict attention to purely 
monotone schemes, then only linear interpolation can be used. 
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However, in [6], it is noted that the solution of the PDE (1.3) can be approximated by the 
solution of a switching system of PDEs with a small switching cost. There, it is shown that the 
solution of the switching system converges to the solution of (1.3) as the switching cost tends to 
zero. In [6j, the switching system was used as a theoretical tool to obtain error estimates. 

Building on the work in |6], the main results of this paper are: 


We formulate the PCPT algorithm in terms of the equivalent switching system, in contrast 
to 127!. We then show that a non-monotone interpolation operation between the switching 
system meshes is convergent to the viscosity solution of (1.3). The only requirement is that 
the interpolation operator be of positive coefficient type. This permits use of limited high 
order interpolation or monotonicity preserving (but not monotone) schemes. 


• We will include two numerical examples. The first example is an uncertain volatility prob¬ 
lem [301 [2] with a bang-bang control, where we demonstrate the effectiveness of a higher 
order (not monotone) interpolation scheme. The second example is a continuous time mean- 
variance asset allocation problem m- In this case, it is difficult to determine the optimal 
policy at each node using analytic methods, and we follow the usual program of discretiz¬ 
ing the control and determining the optimal value by exhaustive search. We compare the 
numerical results obtained using PCPT and a standard policy iteration method. Compara¬ 
ble accuracy is obtained for both techniques, with the PCPT method having a considerably 
smaller computational complexity. 


1.2 Outline 


In order to aid the reader, we provide here an overview of the steps we will follow to carry out our 
analysis. We will write the PCPT algorithm in the unconventional form 


V, 


71+1 


— max 


k=l 




Ar 


- = 0 , 


(1.5) 


where the optimization step is carried out at the beginning of the new timestep, as opposed to 
the conventional form whereby the optimization is performed at the end of the old timestep. Note 
that the scheme is a time-implicit discretization for each fixed control qj , while the optimization is 
carried out explicitly. As discussed above, a decided advantage of this approach is that this scheme 
is unconditionally stable and yet no nonlinear iterations are needed in every timestep. 

In order to carry out our analysis, we perform the following string of approximations: 



HJB equation 



Control discretization 




Switching system 




Discretization 


Vr — sup LgV = 0 
q&Q 

Vr — max = 0 

Qj&Qh 

- Lq.Vj, Vj - i:nax(I4 - c)) = 0 





( 1 . 6 ) 

(1.7) 

( 1 . 8 ) 


(1.9) 
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In the HJB equation (1.6), the control parameter q is assumed to take values in a compact set Q, 
and for fixed q, Lg is a second order elliptic operator as per (1.2). The compact set is discretized 
by a finite set Qh = {?i,. .. j 9j}, where H is the maximum distance between any element in Q 
and Qh- Of course, in the case of a bang-bang control, the admissible set is already discrete. 

The resulting equation (1.7) can then be approximated by the switching system (1.8) as in [6] 
when the cost c > 0 of switching between controls j = 1,..., J goes to zero. When c —)> 0, every Vj 
converges to the solution of (1.7). The switching cost is included to guarantee that the system (1.8) 


satisfies the no-loop condition [2Tj, and hence a comparison property holds. We then freeze the 
policies over time intervals of length At, i.e., restrict the allowable policies to those that assume 
one of the qj over such time intervals, and discretize the PDEs in space and time. Here, we use the 
same timestep Ar for the, say, backward Euler time discretization of the PDE, but generalizations 
are straightforward. We provide for the possibility that the PDEs for different controls are solved 
on different meshes, and in that case interpolation of the discretized value function 14 for control 
qk onto the j-th mesh is needed. We denote by 14 q) this interpolant of 14 evaluated on the j-th 
mesh. 

The remainder of this article is organized as follows. We conclude this section by giving standard 
definitions and assumptions on the equation (1.1). Section shows that the control space can be 
approximated by a finite set, which prepares the formulation as a switching system. Section 
introduces a discretization based on piecewise constant policy timestepping, while Sectionj^contains 
the main result proving convergence of these approximation schemes satisfying a standard set of 
conditions, to the viscosity solution of a switching system. Section [^constructs numerical examples 
for the mean-variance asset allocation problem and the uncertain volatility option pricing model. 
Section [6] concludes. 


1.3 Preliminaries 

We now give the standard definition of a viscosity solution before making assumptions on F. Given 
a function / : D — )• M, where H C M"" open, we first define the upper semi-continuous envelope as 

/*(x) = lim sup{/(y) | yGH(x,r)nH}, (1.10) 

where H(x, r) = {y € | |x — y| < r}. We also have the obvious definition for a lower 

semi-continuous envelope /*(x). 

Definition 1 (Viscosity Solution). A locally bounded function 1/ : H —)• M is o viscosity subsolution 
(respectively supersolution) of 0 if and only if for all smooth test functions 4> £ , and for all 

maximum (respectively minimum) points x of U* — 4> (respectively — (j)), one has 

F* (x, H*(x),D</>(x),D^(/)(x),f7*(x)) < 0 

^respectively F* (x, [/*(x), D i?i(x), (/)(x), [/*(x)) > 0^ . (1-11) 

A locally bounded function U is a viscosity solution if it is both a viscosity subsolution and a viscosity 
supersolution. 

Remark 1 (Smoothness of test functions). Definition^specifies that (f> G C°°, whereas the common 
definition uses cj) € C'^. The equivalence of these two definitions is discussed in Letting 

cj) G C°° simplifies the consistency analysis. 
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Assumption 1 (Properties of F{-)). We assume that Q is a compact set and that aqa"^, fj,q,rq, fq 
are bounded on x Q, Lipschitz in x uniformly in q (i.e., there is a Lipschitz constant which 
holds for all q) and continuous in q. 

Remark 2 (Comparison principle). Assumption^is the same as the one made in It guarantees 
(see, e.g., mm) that a strong comparison principle holds for F, i.e., if V and W are viscosity 
sub- and super solutions, respectively, of El, with V{-,0) < then V <W everywhere. It 

also ensures the well-posedness of the switching system ( ff.4 ), see l2^ and m- 


2 Approximation by finite control set 


In this section, we analyze the validity of the first approximation step from (1.6) to (1.7), i.e., that 
the compact control set may be approximated by a finite set. More precisely, for compact Q C M™' 
(i.e., m G N is the dimension of the parameter space) and Qh C Q such that 


max min \q — q\ < h, 
q&Q q&Qh 


we define a discrete control HJB equation by 


0 = Fh{x,V,DV,D^V) 


Vr - SUPggQ^ LqV, 

v(x)-g(x), 


X X (0,T], 

X G X {0}. 


( 2 . 1 ) 


( 2 . 2 ) 


The following lemma will be useful. 

Lemma 1 (Properties of F(-)). Under Assumption^ for any x and any C°° test function cj) 

|T(x, 4>{x),D4){x),D‘^(j){x)) - Fh{x, (j){x) + D(j){x), D‘^(j){x))\ < w/,(x, h) + 

where 0 Jh{x, h) —)• 0 as h —)• 0, 

w^(^) —)• 0 as ^ 0, (2.3) 


where 0 Jh{x,h) is locally Lipshitz continuous in x, independent of h. 

Proof. See Appendix [A| □ 

Theorem 1. Given Assumption\^ let V and 14 be the unique viscosity solutions to 

F{x,V,DV,D^V) = 0, (2.4) 

Fh{x,Vh,DVh,D^Vh) = 0. (2.5) 

Then I 4 —)• P uniformly on compact sets as h —>■ 0. 

Proof. A consequence of Lemma is that 

•Pfe(x,'^(x)+ ^,...) < F*{x,(j){x),...)+uJh{x,h)+oj^{f,). (2.6) 


Define, for all x. 


Zh(x) 


liminf Vfe(y), 


(2.7) 


Z(x) 


liminf P^(y). 
h^O 
y-^-x 


( 2 . 8 ) 
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We claim that y is a viscosity supersolution of equation (2.4). To show this, fix x and choose 
a smooth test function cj) such that x is a global minimum of y — (/>, and that 

Z(x) = (/>(x). (2.9) 

Then, there exists a sequence x/j —)• x, /i —)> 0, ]d^(xh) )K(x), such that is a global 

minimum of V^f^{xh) — 4>{^h)- At each point x/i, since 14 is a viscosity solution of (2.5), 

0 < Fj^{±h,Vh(.^h),D4>{^h),D‘^4'{^h))- (2.10) 


Let YLhi^h) = 4’i'^h) + Ch, ■Ch 0, /i —>■ 0, so that equation (2.10) becomes 

0 < F^{xh,4>{^h) + ^h,D(j){xh),D^4>{xh))- 
From equations (2.6) and ( |2.11 ) we obtain 

0 < F*{xh,(j){xh),D4){xh),D‘^4i{xh))+uJhi:>^,h)+uj^{^h), 


( 2 . 11 ) 


( 2 . 12 ) 


which gives us 


0 < limsup F*{xh, 4>{^h), D(p{xh), D‘^4>{xh)) + limsup (a;fe(x,/i) + a;^(4)) 

h^O ; 

Xh—X—>x 

< F*{x,4>{x),D(j){x),D‘^4>{x)) 

= F*{x,V{x),D(j){x),D^(l){x)), 


(2.13) 


where we use equation (2.9) and Lemma 


Using similar steps, we can show that V defined similar to (2.8) is a viscosity snbsolution of 


equation (2.4). Invoking the strong comparison principle, V = U = U. Uniform convergence on 
compact sets follows using the same argument as in Remark 6.4, of |16j . □ 


3 Piecewise constant policy timestepping 

Here, we use the equivalence of (O) and ( |1.8[ ) established in [6] as c 
precisely. Consider the HJB eqnation 


0 to formulate (1.9) 


14 = max Lq V, 

qj&Q 

Q = {^1,92, • ■ ■ ,gj}, (3.1) 


where we assume a discrete set of controls Q. We have shown in Section that the optimal value 
under controls chosen from a compact set can be approximated by a control problem with a finite 
set. 


According to [6], we can also approximate equation (3.1) by a switching system. Let Uj,j = 
1,... ,J, be the solution of a system of HJB equations, with 


min 



LqMj, Uj - (max {Uk - c)) 

Uj - g{x) 


0, X E X (0,r], 

0, X E X {0}. 


(3.2) 


7 
















The constant c > 0 is required in order to add some small transaction cost to switching from j ^ k. 
This cost term also ensures that only a finite number of switches can occur (otherwise there would 
be an infinite transaction cost; see also Remark]^. It is shown in [6] that Uj ^ V as c0, for all 

j- 

For computational purposes, we define a finite computational domain C Let 50 denote 
the portions of 50 where we apply approximate Dirichlet conditions. We use the usual notation 
for representing (3.2). Define 

X = (5, r) , DU = {Ur, Us) , D^U = Uss , (3.3) 

and let Bj{x.) be the approximate Dirichlet boundary conditions on 50. Then, the focafeed problem 
is defined as 


0 = Bj {^^B'j^DUj,D‘^Uj,{Uk}kf!=j) 

Uj^r - Lq.Uj, Uj - (maxfc^j {Uk - c)) 


(3.4) 


= 


mm 


X e o\50 X (o,r], 

Dj(x)-^(x), xG 0 X {0}, 

Uj{'x) —Bj{:x), X G 50 x (0,T], 


for j = 1,..., J. Letting pj = DUj, Sj = Uj, we can write equation (3.4) as 

Fj (x, Uj,pj,Sj, {Uk}k^j) = 0. 


(3.5) 


Note that system (3.4) is quasi-monotone (see 123!) , since 

Fj {^,Uj,pj,Sj,{Uk}k^j) < Fj {ic,Uj,pj,Sj,{Wk}k^j) Uk>Wk-, j. (3.6) 


We include here the definition of a viscosity solution for systems of PDFs of the form (3.5) as used 

in [sainiEi]. 

Definition 2 (Viscosity solution of switching system). A locally hounded function [7 : 0 —)• is a 

viscosity subsolution (respectively supers olution) of (3.5) if and only if for all smooth test functions 
4>j G , and for all maximum (respectively minimum) points x ofU* — (fj (respectively Uj^ — 4>j), 
one has 


Fj^ (X, 


LT(x),D</>j(x),D^0j(x),{[/fc(x)}fc/j) < 0 

respectively F* (x, Uj^{x.), D (fj{^), D'^ {C^fc*(x)}fc^j) > 0 | . 


(3.7) 

A locally bounded function U is a viscosity solution if it is both a viscosity subsolution and a viscosity 
supersolution. 

Remark 3. Note that the j-th test function only replaces the derivatives operating on Uj. The 
terms which are a function of Uk,k j, are not affected. 


We discretize (3.4) using the idea of piecewise constant policy timestepping. Define a set of 
nodes Sjj and timesteps r"’, with discretization parameters h and At, i.e.. 


max min 5 — 5,- 4 = h, 

l<j<J i 

Sen 

max(r"+^-r") = Ar. 


(3.8) 











The distinction between At and h is useful for the formulation of the algorithm, but somewhat 
arbitrary for the analysis. We will therefore label meshes and approximations by h and assume 
that 

At —7> 0 as /i —)• 0. 


Define 


= (3.9) 

where r”) refer to points on a specific grid j, and the set of grid points on the grid parame¬ 
terized by h is 

Then we denote the discrete approximation to C/j(x”j) on a grid parameterized by h by Uj{h, x"j), 
which is extended to a function Uj{h, •) on D x {t„} by interpolation. We will sometimes use the 
shorthand notation 


li = Uj {h, xT), xT = {Sj^i , r”), 


(3.10) 


where the dependence on h is understood implicitly. Note that by parameterizing x”j by the control 
index j, we are allowing for different discretization grids for different controls. 


Let £(1. 


be the discrete form of the operator Cq.. We discretize equation (3.4) for x G x 

(0,T], using At = — r” constant for simplicity, by piecewise constant policy timestepping 


u 


n+i 


u 


J,* 


J,* 


where 


= max 


- ArL(l,u?+^ 






— c 


= u 


n-\-i 

3^ 


j = 1,..., j. 


(3.11) 


is the value of this interpolant Uk{h, ■) of Uk{h,x^-) at the z-th point of grid ^j^h- 


(3.12) 


Discretization (3.11) applies the “max” constraint at the beginning of a new timestep. Conven¬ 


tionally, one thinks of piecewise constant policy timestepping as applying the constraint at the end 
of a timestep. 


n+i 


Clearly, these would be algebraically the same thing if Uj ^ ^ instead of uT was considered as 
approximation to Uj{xj-). So at the final timestep, these two possible approximations only differ 
by a final max operation. However, it will be convenient to apply the constraint as in equation 


(3.11). We can then rearrange equation (3.11) to obtain an equation in the form 




n+l 


or b^n / 


= mm 


b^r 

- max - c 


At 






= 0, X G X (0,T]. 


(3.13) 


In the event that Lq is strictly elliptic, then we can interpret the effect of switching at the beginning 
of the timestep as adding a vanishing viscosity term to the switching part of the equation. However, 
note that we do not in general require that Lq be strictly elliptic. 
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We omit the trivial discretizations of Fj{-) in the remaining (boundary) portions of the compu¬ 
tational domain. Note that the notation 


or b^n 


(3.14) 


refers to the set of discrete solution values at nodes neighbouring (in time and space) node {i^n + 1). 

4 Convergence of approximations to the switching system 

Here, we prove convergence of the piecewise constant policy approximation to the solution of 
the switching system (1.8). We start by summarizing the main conditions we need. We will verify 


that these conditions are satisfied for all the schemes we use in our numerical examples. 

Condition 1 (Positive interpolation). We require that the interpolant of the k-th grid onto 

the i-th point of the j-th grid can be written as 


u 


Q;GW^(_7,2,n) 


^k,iU),aMh,^k,c 


(4.1) 


where 




U 


k,i{j),a 


> 0 , 


(4.2) 


and N^{j,i,n) are the neighbours to the point on grid ^lk,h- 


Remark 4 (Monotone versus limited interpolation). Condition (f-l), (4.2) is elearly satisfied by 


linear interpolation on simplices and multi-linear interpolation on hyper-rectangles, in which case 
the weights ^ki{j)a only functions of the coordinates j and This interpolation is then 

also a monotone operation and results in an overall monotone scheme, as defined below in Condition 

H 

We can, however, relax the requirement of monotonicity of the interpolation step by allowing 
weights ^ki{j)a ~ ^k i{j) a^'^V)> possibly nonlinear functions of the interpolated nodal values. 
It is easy to see that as long as 


min, < K,i(j) < max 


(4.3) 


one can then always find weights such that (f.D, {4-^) hold. One can enforce (4-3) by o, simple 


limiting step applied to any, potentially higher order, interpolant. 


Generally, the decomposition (4.1) will not be the way by which a particular interpolation is 


originally defined or constructed in practice. The point is that whenever (4.3) holds, weights of this 
form exist and this is all that is needed for the theoretical analysis. An example are one-dimensional 
monotonicity preserving interpolants such as those in m- As discussed in m, the interpolation 
in 121]/ is eonstructed to be monotonicity preserving (i.e., the interpolant is increasing over intervals 


where the interpolated values are increasing), and therefore satisfies condition (4.2) 


10 











Condition 2 (Weak Monotonicity). We require that discretization (3.13) he monotone with respect 


to i.e., if 



> 


A{i,j,n), 



> 


\/{i,k,n), 

(4.4) 

then 

\ or b^n 

fc/j 

) <G, | 

\ orb^n / 

1 (4.5) 


Remark 5. Note that the requirement that the scheme he monotone in u^- (the interpolated solu¬ 
tion) is a weaker condition than requiring monotonicity in In particular, let us re-iterate that 
we are not requiring interpolation to he a monotone operation, as long as Condition^is satisfied. 


Condition 3 {loo stability). We require that the solution of equation (3.13), Uj{h,x^f^), exists 
and is bounded independent of h. 

Remark 6. The bounds (4-3) implied by Condition^ ensure that the interpolation step does not 
increase the loo norm of the solution. 

Condition 4 (Consistency). We require that we have local consistency, in the sense that, for any 
smooth function fij, and any functions pk (not necessarily smooth) 


G 


\ or b^n / 


< UJi{h) +UJ2{0, 

u}i{h) —>-0 as /i —>■ 0, 0 J 2 {fi) —)• 0 as ^ 0. (4.6) 

Theorem 2. Under Assumption^ the solution of any scheme of the form (3.13) satisfying Con¬ 
ditions M converges to the viscosity solution of {3.4% uniformly on hounded domains. 

Proof. We basically follow along the lines in [7] , with the generalizations in m for weakly coupled 
systems. We include the details here in order to show that we only require Condition for the 
interpolation operator, and this permits use of certain classes of high order interpolation. 

Dehne the upper semi-continuous function u by 


Uj{x.) = limsup Uj{h,x^f^) , 
h^o 

n + 1 


(4.7) 






where G ^j,h- Similarly, we define the lower semi-continuous function u by 

uAyi) = liminf Uj(/i, . 


h-s.0 

..n+i , 


(4.8) 


Note that the above definitions imply that u* = Uj and = u^. 
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Let X be fixed and (pj be a smooth test function such that 

4>j{yi) = Uj{x), 

(pji'x.) > Uj{x.), X X. (4.9) 

This of course means that {uj — 4>j) has a global maximum at x = x. Consider a sequence of grids 
with discretization parameter hi, such that —)> 0 for I —)• oo. We use the notation 

E (4.10) 

to refer to the grid point + 1) on the grid parameterized by hi, associated with discrete 

control Qj. Let be the point on grid ^j^hi such that 




(4.11) 


has a global maximum, where (f>j is a test function satisfying equation (4.9). Note that in general, 
for any finite /i£, ^ x. 

Following the usual arguments from [7], and more particularly in for a sequence of grids 


^j,he parameterized by hi, there exists a set of grid nodes {ie,ni + 1), such that (4.11) is a global 
maximum and 


x?^+^ X, xjj^ X, Uj{he, x”j+^) Uj{x), for oo, 


J,ie 


where E ^j,he and, for k ^ j, noting equation (4.2) for the interpolant 

have 


u 


ni _ 


k,ieij) 




^k,ii(j),ai^k,ai ~ 


ofeAf'=(-) 

hmsup'Ufc^j^(j) < Ufc(x) 
i^OO 


where x^ ^ ^k,hi and x = (5,r). Let 


u 


JM 




4f‘ - 


JM 


k^j , 


oo , 
(4.12) 


(4.13) 


Note that for any finite mesh size hi the global maximum in equation (4.11) is not necessarily zero, 
hence we define by 


ne+l 




such that 


6 > 0 


for i 


oo . 


(4.14) 
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then 


Since (4.11) is a global maximum at 


} 0.1^11 < ■ 
or hi^ni or hii^ni 




(4.15) 


Substituting equations (4.14) and (4.15) into equation (3.13), and using the monotonicity prop¬ 
erty of the discretization (4.4) gives 


0>Gj{ 


or bj^n 


(4.16) 


Note that we do not replace by the test function in equation (4.16). This is because the 

test function is only defined such that ipj > Uj, and there is no such relationship with Uk,k ^ j. 
Let 


Pk% = Pki^j'ie^ = max({ifc(x^.^),nfc(i)) . 


From equation (4.12), we have that 


and that 




PkM - 


(4.17) 

(4.18) 

(4.19) 


discretization (4.4) gives 


Substituting equation (|4.17|) into equation (|4.16|), and using the monotonicity property of the 

(4.20) 


o>G4x”jy 


. 'x. {<*52') •#. + {pZ,}i=*i) 

or b^n / 


Equations (|4.20|) and (|4.6|) then imply that 

n£+l\ 7-,2 


0 > Fj{ X 


nt+l 


> ^), D (()j(x”^"), ^j(x”"T"), {pk{x]^)}k^j 1 - uji{he) - u}2{^e) . 


nt+l 


Recalling equation (4.18), we have that 


0 > liminf Fj ((>j(x"/+^), D (/>j(x”/+^), (^j(x"/+^), 


— lim inf uJi{hi) — lim inf W 2 {£,{) 

i^OO i^OO 


> Fj^ y^,4>j (x), T) </>j (x), 4>j (x), {Uk (x)} 

= Fj^(^,Uj{x),D(pj{x),D‘^ 0j(x), {uk{x)}kjLj^ 


(4.21) 


Hence uj is a subsolution of equation (3.4). A similar argument shows that Uj is a supersolution of 
equation (3.4). Since a strong comparison principle holds for the switching system under Assump¬ 
tion (see Proposition 2.1 in [6] and Remark [^, we then have Uj = Uj is the unique continuous 
viscosity solution of equation ( |3.4[ ). 

We have thus shown the result. □ 
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Remark 7. 

- K,,iyj)- -- “ - _|- 

a problem in terms of consistency, since we cannot assume that w 


Note that we have j(j)} appearing in equation (3.13), which would appear to cause 


n 

k,a 


111 

are not necessarily smooth. The key fact here is that since equation \ 4 .I 8 ) holds, we do not need 
smoothness. 

Remark 8 (Interpolation and switching cost). We include here a brief example of the role of the 


switching cost c in (3.2). Assume we solve the degenerate equation ut = 0 and write it as the trivial 
HJB equation supQgjQ = 0. We represent and vf on two different meshes with nodes at 
xf = (i + Q.h9)h, i.e., shifted by half a mesh size. The fully implicit discretization over a single 
timestep is = u®’"" for both controls. Consider now the situation of convex {u^), such that 

linear interpolation increases the solution. Using piecewise linear interpolation with c = 0 at the 
end of the timestep. 


u- 


( u 

= 


+ u. 


i+l 


^)/2, 


and similarly for . Repeating this for the next timestep, the piecewise constant policy dis¬ 

cretization is equivalent to 


d,n+2 1 

Uf =|u, 


6,n 

i-l 


, 1 6,n . 1 O.n 

+ 2“i +4^4 i 




6 , 71+2 9,71 

+ - + 


6,71 

^i -1 


6,71 

h 


6,71 

h+1 


ih+y 


/l2 


If we pick At = hf /S, this discretization is consistent with the standard heat equation ut = Ux 
instead ofut = 0. 

With the cost c > 0 switched on, for sufficiently small h, we will have 


u. 


l.n , / 0,n+l , 0,n+l\ /„ 

>{uf )/2-c, 


such that and, by induction, the equation ut = 0 is solved exactly. 

For an overall convergent method, one has to pick h and At as a function of c, depending on 
the interpolation method and smoothness of the solution, and then let c ^ 0. We discuss this in 
more detail in Section\^for concrete examples. 


Remark 9 (Reduction to standard piecewise constant policy method [271 [6]). Discretization (3.11), 
with c = 0, can be viewed as a form of the usual piecewise constant policy method As a 

result, if linear interpolation is used to transfer information between grids, then (3.11) is a monotone 


discretization of HJB equation (3.1), which is easily shown to satisfy the standard requirements for 
convergence to the viscosity solution. If we use standard finite difference schemes, we would expect 
the spatial error (for smooth test functions) to be of size 0{h) with timestepping error of size 0(At). 
In this case (c = 0), if the solution is smooth, we would expect the total discretization error to be 
of size 0{h) + 0(At) + Ofhf / At), where the last term arises from a linear interpolation error 
accumulated 0(1/Ar) times. Note that this term will be absent in the case c > 0, since the finite 
switching cost will prevent an 0(1/At) accumulation of interpolation error. 


5 Numerical examples 

In this section, we study the convergence of discretization schemes based on piecewise constant 
policy timestepping in numerical experiments. We present two examples, the uncertain volatility 
model from derivative pricing, and a mean-variance asset allocation problem. In both examples, 
we investigate the convergence with respect to the timestep and mesh size. 
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r 

min 

^ max 

T 

K 

Ki 

K2 

So 

0.05 

0.3 

0.5 

1 

100 

80 

120 

100 


Table 5.1: Model parameters used in numerical experiments for uncertain volatility model. 


5.1 The uncertain volatility model 

We study first the uncertain volatility option pricing model [30]. In this example, we also examine 
the role of the switching cost and the impact of different interpolation methods. 

The super-replication value of a European-style derivative is given by the HJB equation 



dV 

- -sup L^V = 0, 

OT o-es 

(5.1) 

where 

1 n n(92l/ dV 

L.V = \a^S^^+rS—-rV 

(5.2) 

for S G (0, oo) 

, r G (0, T], with t = T — t and 



— [^min;'^max] • 

(5.3) 


In addition to the PDE, the value function satisfies a terminal condition. Eor the numerical tests, 
we choose a butterfly payoff function P such that 

V{S, 0) = P{S) = max(S' — it'i, 0) — 2 max(5 — K,0) max(5 — K2,0) 

and we localize the domain to [5min, S'max] with 

V{Sra..,T) = 0, (5,r) G {5„iax} X (0, T] , 

w = -rV, (5,r) G{5a,in} X (0,r], 


and parameters as in Table 5.1 

It is well understood m that the optimal control is always attained at one of the interval 
boundaries (‘bang-bang’), depending on the sign of the second derivative of the value function V. 
The payoff function was chosen such that V has mixed convexity and therefore the optimal control 
differs between different regions of the state space and changes in time. 


We use logarithmic coordinates X = logS, so instead of (5.2) we approximate 


L^V = 




a 




(5.4) 


and now the coefficients are bounded and X G [log 5min, log S'max]) so Assumption is satisfied. 
Here, it is straightforward to construct “positive coefficient” schemes for the linear PDE arising 
when the control set is a singleton. A well-established route allows us to construct monotone, con¬ 
sistent and loo stable schemes for the fully non-linear problem from positive coefficient discretization 
operators using a direct control method of the “discretize, then optimize” type, see [20] . 


/ yn+l _ yn \ 

min f---=0, (5.5) 

ie{i,2} V Ar y ’ 


15 
















where qi = cJmin and q 2 = (Tmax- 

The resulting non-linear finite dimensional system of equations can be solved by policy iteration 
(Howard’s algorithm, see [9]). In particular, we will use standard central finite differences in space 
and an implicit Euler discretization in time. For small enough h, this scheme is monotone, ioo 
stable and consistent in the viscosity sense. 

We compare the direct control method to several variants of piecewise constant policy timestep¬ 
ping on the basis of these discretizations, as described in Section]^ specifically equation (3.11). 
Here, only two control values have to be considered and the switching system is two-dimensional. 
In the case of no interpolation, the scheme simplifies to 


- maxfcg{i^2} {Vk - <^k,j) 

At 


rh -fT-n+l 

’ 


Ck,j — 


j = l,2. 



(5.6) 


We again discretize Lq. using a positive coefficient discretization, hence it is straightforward to 
verify that Assumptionand Conditions [T||^ hold (on the localized domain S G [Smin, *S'max])- 
A solution extrapolated from the finest meshes was computed as an approximation to the exact 
solution and used to estimate the errors. The numerical value of this solution is V{So, 0) = 1.67012 
(see also [33]). From this, we approximate the error as 

e(/i, Ar) = |H(5'o,0) - t/i(/i, At, c; 5o, 0)|, (5.7) 

where V is the exact solution and [/i(/i, Ar, c; •, •) the numerical approximation to Ui, the first 
component of the switching system, for mesh size h, timestep Ar, and switching cost c. 


Dependence on timestep At and switching cost c 


We first analyze how the switching cost affects convergence of the approximations. In Fig. 5.1 
compare the following two cases. 


we 


1. Policy timestepping, fixed mesh: We use (5.6) on a single uniform mesh on [log(Ar) — 
4 • cj,log(it') + 4: ■ a], i.e., encompassing four standard deviations either side, where a is the 
average of the two extreme volatilities. 

For fixed switching cost, the error first decreases linearly as the timestep decreases, but 
eventually converges to a non-zero value. This convergence is monotone. Moreover, for fixed 
(sufficiently small) h, the asymptotic error for At —)• 0 is decreasing with c. 

2. Policy timestepping, linear interpolation: We also study the use of separate meshes for 
the two components of the switching system, Vi and V 2 - In particular, we use uniform meshes 
on the intervals [log(A:) - 4 • amin, log(Ar) + 4 • CTmin] and [log(A:) - 4 • Umax, log(A') -h 4 • Umax], 
so that for the chosen parameters mesh points on the two meshes do not coincide. Then, 
interpolation is necessary to represent these solutions on both meshes and to evaluate the 
explicit terms (those at time-level n -|- 1) in (5.6). In the case of linear interpolation, the 
overall scheme is monotone, and trivially satisfies Condition]^ As a result of the switching 
cost, the cumulative effect of linear interpolation in each timestep is controlled even for small 
At (see Remark]^. 
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Varying Cost, No Interpolation 


Varying Cost, Linear Interpolation 




Figure 5.1: 
for different methods the log 2 error, 


timestep (At)^ = 1/8-2 


The uncertain volatility test case with parameters as in Table 5.1 


Shown is 


where Cm = e((Ar)m,/i, c) from (5.7[ ) is the error for 
{1/8,..., 1/32768} and in each plot, from top to bottom. 


c = 1/20,1/40,1/80,1/160,1/320,1/640. The mesh size is fixed ath = 1/1024. The dotted line has 
slope —1. Left: Piecewise constant policy timestepping on a single mesh; right: linear interpolation 
between individual meshes for each control. 


Next, we analyze the convergence jointly in h and At for fixed c, as well as the convergence 
with respect to c for the case with interpolation. The results are given in Table 5.2 For fixed 


positive switching cost c > 0, we compute approximations on a sequence of time and space meshes 
with iVfc = 2iVfc_i and = 2Mk_i. The asymptotic ratio of about two is consistent with an error 
of 0(h) + 0 (At). 

We now study the difference between the solution of the switching system for fixed cost c 


and the solution of the HJB equation (5.1). Considering the boldface values in the table as good 
approximations for this difference, we observe convergence as c —)• 0 which is roughly consistent 
with order 3/4. The theoretically proven order of 1/3 from Theorem 2.3 in [6] does not seem sharp 
for these data. 

The approximation error in At and h does not appear to be affected by c as long as c > 0, 


which is seen by comparison of the lines ‘(c)’ in Table 5.2 for c = 1/10 to c = 1/640, for small 
enough mesh parameters (last three columns). In computations, it would therefore seem prudent 
to pick At = 0{h) and c = i.e., proportional errors, even though faster, but more erratic 

convergence is obtained setting c = 0 uniformly. 


Dependence on timestep At and mesh size h 


In Fig. 5.2, we show the convergence in both the timestep and mesh size for two different costs, 
c = 0.01 and c = 0.16. Compare this to the case c = 0 in Fig. 5.3 (middle row, left). We can 


make a number of observations. For large switching cost (right plot), the difference between the 
switching system and the HJB equation is large and dominates the discretization error. For fixed c 
and h, there appears to be convergence as At —0, and if we also let h —0 the solutions converge 
to the solution of the switching system with fixed c. For comparable values for h and At, there 
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Nk 

Mk 


32 

512 

64 

1024 

128 

2048 

256 

4096 

512 

8192 

1024 

16384 

2048 

32768 

4096 

65536 

c 










1/10 

(a) 

2.0692 

1.9724 

1.9660 

1.9532 

1.9491 

1.9478 

1.9474 

1.9472 


(b) 

0.3991 

0.3023 

0.2958 

0.2831 

0.2789 

0.2777 

0.2773 

0.2771 


(c) 


-0.0968 

-0.0065 

-0.0127 

-0.0042 

-0.0012 

-0.0004 

-0.0002 


(d) 



14.9669 

0.5079 

3.0621 

3.4257 

2.8297 

2.3076 

1/40 

(a) 

1.6126 

1.5441 

1.7406 

1.7663 

1.7623 

1.7612 

1.7608 

1.7606 


(b) 

-0.0575 

-0.1260 

0.0705 

0.0961 

0.0922 

0.0910 

0.0906 

0.0905 


(c) 


-0.0685 

0.1965 

0.0256 

-0.0040 

-0.0011 

-0.0004 

-0.0002 


(d) 



-0.3486 

7.6692 

-6.4629 

3.5076 

2.7965 

2.3037 

1/160 

(a) 

1.1989 

1.2205 

1.5510 

1.7019 

1.7033 

1.7022 

1.7018 

1.7016 


(b) 

-0.4712 

-0.4496 

-0.1191 

0.0318 

0.0331 

0.0320 

0.0317 

0.0315 


(c) 


0.0216 

0.3305 

0.1509 

0.0013 

-0.0011 

-0.0004 

-0.0002 


(d) 



0.0653 

2.1902 

112.7836 

-1.2152 

2.7969 

2.2994 

1/640 

(a) 

0.9256 

0.9897 

1.3752 

1.6448 

1.6833 

1.6822 

1.6818 

1.6816 


(b) 

-0.7446 

-0.6804 

-0.2949 

-0.0254 

0.0132 

0.0121 

0.0117 

0.0115 


(c) 


0.0641 

0.3855 

0.2696 

0.0385 

-0.0011 

-0.0004 

-0.0002 


(d) 



0.1664 

1.4301 

6.9932 

-35.1108 

2.7835 

2.3009 

0 

(a) 

0.8271 

0.7221 

1.0227 

1.4109 

1.6654 

1.6702 

1.6703 

1.6702 


(b) 

-0.8430 

-0.9480 

-0.6474 

-0.2592 

-0.0047 

0.0001 

0.0002 

0.0001 


(c) 


-0.1050 

0.3005 

0.3882 

0.2545 

0.0048 

0.0001 

-0.0001 


(d) 



-0.3494 

0.7742 

1.5252 

53.3568 

39.0446 

-1.6580 


Table 5.2: The uncertain volatility test case with parameters as in Table 5.1: piecewise constant 


policy timestepping with linear interpolation between individual meshes for each control; convergence 
with respect to mesh parameters and switching cost. Shown are: (a) the numerical solution Vk = 
V{Nk, Mk,c; So,0); (b) the difference to the exact solution Vk — V{So,0); (c) the increments Vk — 
Vk-i; (d) the ratios of increments {Vk — Vk-i)/{Vk-i — Vk-2)- 
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Fixed Cost, Mesh Refinement, Linear Interpolation Fixed Cost, Mesh Refinement, Linear Interpolation 




ID 12 14 16 


Figure 5.2: The uncertain volatility test case with parameters as in Table 5.1 Piecewise 
constant policy timestepping with linear interpolation between individual meshes for each con¬ 
trol. Shown is the log 2 error, where Cm = 

(Ar)^ = 1/8 • 2-™ e {1/8,..., 1/524288} and 


in 


h,c) from {5.1) is the error for timestep 
each plot, from top to bottom, h = 


1/4,1/16,1/64,1/256,1/1024,1/4096. The cost size is fixed at c = 0.01 (left) and c = 0.16 (right). 
The dotted line has slope —1. The downward spikes are a result of error cancellation for a particular 
combination of h, Ar and c. 


appears to be a cancellation of leading order errors in h and At with opposite signs, which appears 
as downward spikes in the left and middle plot. 

In this particular case of linear interpolation, since the overall scheme is monotone, convergence 


to the viscosity solntion of (5.1) is ensnred if c = 0 as long as the discretization is consistent. Recall 
from Remark]^ that (for smooth solutions) the discretization error is of the form 0(Ar) + 0{h?) 
+ 0{h‘^/A t), with the third term being the cumulative effect of linear interpolation over 0(1/At) 
timesteps. Hence we can ensure consistency by requiring that At = 0{h). 

We now analyze in more detail the convergence in At and h for the degenerate case with c = 0. 


The computational results are shown in Fig. 5.3 and are discnssed in the following list 


1. Policy timestepping, fixed mesh: We use (5.6) on a single uniform mesh on [log(iF) — 
4 • cj,log(iF) + 4 • O'], i.e., encompassing four standard deviations either side, where a is the 
average of the two extreme volatilities. For fixed mesh size, the error approaches a constant 
level for decreasing time-step, but for simultaneously diminishing mesh size the observed time 
discretization error is clearly of first order in the timestep. 

2. Direct control: Here, the optimal control is implicitly found with the solntion as described 


above - see, in particnlar, (5.5) - and therefore we can disentangle the Euler discretization 


error from the effect of piecewise constant control. Comparing the envelope to the curves, 
parallel to the dotted line with slope minus one, shows first order convergence as in the 
previous case, but with a lower intercept which indicates that the time discretization error 
is about a factor 4 smaller. Although the number of policy iterations per time-step was 
consistently small (usually 2-4), the computational time here was dramatically larger (due 
to the need to generate new matrices in each iteration) and therefore solutions could not be 
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computed for the same number of timesteps as for the other cases. 


3. Policy timestepping, linear interpolation: Again, first order convergence in the timestep 
is observed, however, the leading error terms are, from Remark]^ of the form 0(/i) + 0(Ar) + 
0{h?‘/ At). As a result, convergence is ensured only if h goes to zero faster than t/t, and 
for h ~ Ar first order convergence is expected. Because of the maximum norm stability and 
linear interpolation, the error does not explode even as At —)• 0 for fixed h. In fact, the 
solution goes to zero here as the interpolation introduces increasing artificial diffusion (see 
also Remark]^ and the solution is absorbed at the boundaries. 

4. Policy timestepping, linear interpolation, reference mesh: Although not an issue 
for J = 2 control parameters, if the dimension of the switching system is J, the number of 
interpolations from each mesh onto all other meshes is an 0{J^) operation. We can avoid 
this using a single ‘reference mesh’ to keep track of the solution. So in addition to the two 
meshes associated with Vi and V 2 as above under item[^, we introduce a reference mesh, 
uniform on [log(A') — 4- a, log(Ar) + 4 • u], and Vk^i is constructed by linear interpolation from 
14 onto the reference mesh and then onto the i-th mesh point of the j-th mesh. With now 
two interpolations for each solution every timestep, convergence is still of first order but with 
a significantly higher factor than with direct interpolation between meshes. The number of 
interpolations needed for a J-dimensional switching system is now 0{J). 


5. Policy timestepping, cubic interpolation: Finally, to reduce the accumulated interpola¬ 
tion error, we use the possibility of limited higher order interpolation for the mesh transfer 
afforded to us by Condition first without the use of a reference mesh. In particular, we use 
monotone piecewise cubic Hermite interpolation as in [21] . Note that the interpolation |21) is 
monotonicity preserving but not monotone in the viscosity sense [?]• However, Condition 
is satisfied. Note that we use c = 0 in this test case, which, strictly speaking, does not ensure 
convergence to the viscosity solution. However, it is clear from Figure 5.1 that the limiting 
case of c —)• 0 does in fact converge to the viscosity solution, hence it is interesting to include 
the case c = 0. 


The approximation order of the interpolation method in m is guaranteed to be cubic only 
if the data are in fact monotone, and this is not the case for our initial data. Nonetheless, 
the error is significantly reduced compared to the linear interpolation case. 


6 . Policy timestepping, cubic interpolation, reference mesh: The results with cubic 
interpolation onto a reference mesh are not as accurate as for direct cubic interpolation 
between the computational meshes, and have a similar accuracy to the results for linear 
interpolation without reference mesh. 


5.2 Mean-variance asset allocation 

As a second example we study the mean-variance asset allocation problem as discussed in |39j . 
following the embedding technique introduced in |291141j . In this example, we use the same grids 
for each constant policy mesh, and focus on the effects of discretization of the control. This example 
demonstrates that piecewise constant policy timestepping does not introduce any significant extra 
error compared to hrst order Euler timestepping with either known optimal control, or a numerical 
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Policy Timestepping 
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0 



m 


Direct Control 



m 


Linear Interpolation 


Linear Interpolation, Reference Mesh 
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m 


Cubic Interpolation 


Cubic Interpolation, Reference Mesh 
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m 


5.1 


Shown is 
error for 


Figure 5.3: The uncertain volatility test case with parameters as in Table 
for different methods the log 2 error, where Cm = e((AT)m,/i) from (5.1) is the 
timestep {At)^ = 1/8 ■ 2“™ € {1/8,1/524288} and in each plot, from top to bottom, 
h = 1/4,1/16,1/64,1/256,1/1024,1/4096. The dotted line has slope —1. The plots refer to, from 
top left lexicographically: the piecewise constant time-stepping method on a fixed mesh; the direct 
control method; the piecewise constant time-stepping method with: linear interpolation; linear in¬ 
terpolation onto a reference mesh; cubic interpolation; cubic interpolation onto a reference mesh. 
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optimal control obtained implicitly from the finite difference scheme. We will also see that even 
a fairly coarse discretization of the control admissible set yields good results. We have seen this 
property of discretized controls in many examples. 

The method determines the pre-commitment mean variance optimal strategy [ 8 ]. Note that it 
is possible to develop a numerical method for solution of the time-consistent version of this problem 
m- However, since the time consistent problem can be viewed as a constrained solution of the pre¬ 
commitment problem, the time consistent solution is sub-optimal compared to the pre-commitment 
solution. Specifically, we consider here the sub-problem given by the equation 

0, (5.8) 

{w-iy. ( 5 . 9 ) 


dV 

^ - inf CPV = 
OT p£P 


V{W, 0) = 


on (— 00 , 00 ), and ( 0 , 00 ), with 


LPV = laVw‘^ ^ 


d'^V dV 

^ + (7r + »F(r + paf)) — 


(5.10) 


and either p G (— 00 , 00 ) or p G [0,Pmax]- Observe that (5.8) does not satisfy Assumptionsince p 
can be unbounded. We will re-parameterize the control variable to avoid this problem. 


By solving equation (5.8'5.9) for various values of the parameter 7 , we can trace out the efficient 


frontier in the expected value, variance plane 

We use the standard finite difference discretization and make “maximum” use of central differ¬ 
ences |38j whenever a positive coefficient scheme is achieved and use upwind differences only where 
necessary for monotonicity. 


The PDE (5.8-5.9) is specified on an infinite domain. For numerical purposes, we approximate 


this by means of a localized problem, with approximate boundary conditions at finite values of | 1 T|. 
We use an asymptotic approximation of V for large \W\. In the cases we consider, an asymptotic 
value for the optimal control is more directly available than for the value function. Therefore, the 
following solution under constant control will be applied as approximate boundary condition. More 
precisely, the solution to the PDE 




dV 

~d^ 


dV 

mv 


= 0 


(5.11) 


with terminal condition (5.9) is given by 


V(W,t) = a{T)W^ + I3 {t)W+ 6{t), 


(5.12) 


where t = T — t and 

Q((r) = exp((a^-|-26)r), 

/3(r) = -( 7 -I-c) exp(5r)-I-cexp((a^-I-26)r), 

6{t) = + (exp(6T) - 1) -h (exp((o2 -h 26 )t) -1) + '^, where 

0 + 2b 4 

c = 2TT/{a^ + b). 
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r 

a 


TT 

ITo 

T 

7 

A 

0.03 

0.15 

0.33 

0.1 

1 

20 

14.47 

1.762 


Table 5.3: Model parameters used in numerical experiments for mean-variance problem. 


In comparison to [39], who only derive the highest-order term, this gives an asymptotically more 
accurate approximation and allows us to use substantially smaller domains for the computation. 
Following |39|, we use the parameters in Table 5.3 throughout. 


We study two different cases for the permissible sets for state-variable and controls, one where 
VF,p G M, and one where VF > 0, p G [Q,Pmax]- 


Bankruptcy allowed, unbounded control 


If bankruptcy {W < 0) is allowed, the PDF (5 .8 5.9) holds on (— 00 , 00 ). In this case, a closed-form 
solution is known from 


p*(W,t) = - 


, where the optimal policy is given by 

e 


ctW 


W- 




9) 


(5.13) 


Moreover, under this optimal policy, we find from the formulae in 


Var[WT] = 


o-eT 


1 - 


E[Wt\ - ( + 


7r(e^^ — 1) 


E[Wt] = {wo + + 


7(1 




TT 

-e 

r 


-(fT 


such that {^Var[WT],E[WT]) = (0.794,6.784), and E[{Wt - l/2f] = Far[Wr] + E[Wt? - 
'^E\Wt\ + 7^/4 = 0.8338 for the parameters in Table 


5.3 


As the optimal policy in the form (5.13) is unbounded, we perform the control discretization in 
a different control variable. Noting that p*W is bounded as IT —?• 0, it seems natural to consider 
pW as control variable in this area; however, p*W ~ —f^jcrW as |1T| — 00 . This leads us to 
consider 


9 = 


pW 


max(l, w|lT|) 


as control variable for some cj > 0 , and 


~ dV 

C^V = max( 1, VT^) + (^ + ^^ + (?max(l, a;|lT|)cT^)- 


dW 


(5.14) 


(5.15) 


The optimal control q*{W,t) will be bounded on a localized domain, and we fix an interval 
Q = [Qmin,Qmax] in which we search for the optimal control by a crude approximation. In this 
whole process, a precise knowledge of the exact optimal control is not necessary, as we only use the 
rough asymptotic shape. For the computations below, we pick cj = 5 and Q = [—2.5,3.5]. Since 
we solve the PDF on a localized domain (|1T| bounded), and the control is now bounded as well, 
the localized version of equations (5.8-5.9) now satisfies Assumption 
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Figure 5.4: The mean-variance test case with parameters as in Table 5__3 Left: The numerical 
approximation to V (W, 0) for piecewise constant policies for M = 80, TV = 80 and 20 policy 
steps. The dashed line is the asymptotic approximation for \W\ —>■ oo. Right: The corresponding 
approximation to the optimal policy q{W, 0), and the analytical optimal policy {5.13). The horizontal 
lines are at the asymptotic optimal policies for \W\ —>■ oo. 


We note that this is an example where the optimal control is an unbounded function of the 
state variable, but by a suitable reformulation the piecewise constant policy timestepping method 
can still be applied, with the policy chosen from a bounded set. 

From (5.13) one sees that q*{W,t) —)■ —f,/cr for \W\ —)■ oo. Therefore, asymptotically, (5.8), 


(5.10) takes the form (5.11) with suitable a and b, obtained by inserting the constant asymptotic 
optimal policy. We can then use the asymptotically exact boundary conditions (5.12) for both 
Wmax and Wmin- We choose Wmax = 40 and Wmin = —40 in the computations. 

The discretized switching system has the form 


- mini<fc<j( 14 " - Ck,j) 


At 




(5.16) 


where Ckj is defined as in equation (5.6). In this case, we can set the switching parameter Ckj = 0 
since no interpolation is used, and this reduces to conventional piecewise constant policy timestep- 
ping m- Then the numerical approximations to all J components of the switching system are the 
same in each timestep after the minimum is taken. 

Fig. 5.4 shows the value function V and its asymptotic approximation for large |VF|. The two 


functions have visually identical tangents at the boundaries, and indeed experimentation with the 
values of Wmin and Wmax shows that the results around W = 1 are not significantly affected by 
this approximation. 


Also shown in Fig. 5.4 is the approximate optimal policy obtained numerically from the policy 


timestepping discretization with 20 policy steps, and the exact formula (5.13), transformed into a 


bounded control as per (5.14). 


Table 5.4 illustrates the convergence as the control mesh is refined for a fixed time and spatial 
mesh. The estimated order of convergence over these refinement levels is 2. We pick J = 40 fixed 
for the following tests of the convergence in the mesh size and timestep. For this value, the control 
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J = 5 

J = 8 

J = 10 

J= 15 

J = 20 

J = 29 

J = 40 

J = 57 

J = 80 

(a) 

2.257 

1.531 

1.429 

1.254 

1.230 

1.196 

1.186 

1.180 

1.178 

(b) 


-0.725 

-0.101 

-0.175 

-0.0241 

-0.0339 

-0.0104 

-5.4 -10-3 

-2.5 -10-3 

(c) 



7.12 

0.58 

7.26 

0.71 

3.25 

1.92 

2.16 


Table 5.4: The mean-variance test case with parameters as in Table 5.3 Convergence of the 
control discretization alone for N = 480, M = 120. Shown are: (a) the numerical solution T4 = 
V{N,M,Jk;Wo,0); (b) the increments Vk - Vk-i; (c) the ratios {Vk - Vk-i)/{Vk-i - 14-2); for 


Jk — 


5-v^' 


k-l 


Policy Timestepping Exact Policy 




5.3 


Figure 5.5: The mean-variance test case with parameters as in Table |5. Similar to Fig. 
in the previous section, the log 2 error for (At)^ = 1/8 • 2“™ G {1/8, • • ■, 1/32768} and in each 
plot, from top to bottom, h = 1/4,1/16,1/64,1/256,1/1024. The straight line has slope —1. Left: 
Piecewise constant policy timestepping with fO equally spaced policies in [—2.5,3.5]. Right: Using 
the exact policy given by (5.13). 


discretization error was empirically negligible (compared to the time and spatial discretization 
error). 

Fig. |5.5|shows the convergence of the approximations for piecewise constant policy timestepping 


and for the use of the exact policy given by (5.13). In the latter case, the error is solely due to 


the Euler time-discretization and spatial hnite differences. For piecewise constant timestepping, 
40 policies were used, so that the computational time for the same mesh is about a factor of 
40 larger than for a single policy, hence we show slightly fewer rehnements. It appears that the 
spatial approximation error for large mesh size h is smaller if knowledge of the optimal control is 
used. The envelope showing the time discretization error is consistent with hrst order convergence. 
Interestingly, the intercept is about 4 units higher for the exact policy, so that the results using 
policy timestepping are about a factor 16 more accurate for the same timestep. This is not to be 
expected generally and must result from opposite signs of the Euler truncation error and the error 
due to piecewise constant policies. 
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Figure 5.6: 
hounded control. 


The mean-variance test case with parameters as in Table 5.3. with no bankruptcy and 
All model parameters are as in Table 5.3 Left: Value function F(VF, 0) (solid 


line) and the asymptotic approximation (dashed line) for large W. Right: The numerical optimal 
policy p{W, 0). 


No bankruptcy, bounded control 


If bankruptcy {W < 0) is not allowed, the PDE (5.8-5.9) holds on (0,oo). The boundary equation 
at IT = 0 is then 


T-(0,r) - 7rIV(0,r) = 0, 


(5.17) 


see [39] for a discussion. For vr > 0 there is an outgoing characteristic (going backwards in time) 
so that no boundary condition is required and we can approximate ( |5.17 ) by upwind differences 
from interior mesh points. In fact, as we are switching to upwind differences locally whenever the 
monotonicity of the scheme is violated (see above and |38| ), upwinding will always be used for small 
IT if vr > 0. 

For bounded control with no short-selling, P = \f),pmax\ in (5.8). In the computations, we 
choose pmax = 1.5 as an attained upper bound (the same used in [39|). This would correspond to a 


typical leverage constraint. For large IT, we use again the approximation (5.12), with coefficients 


based on the asymptotic optimal control p = 0 (see Fig. 5.6). 


The numerically computed value function (a closed-form solution is not available in this case) 


is shown in Fig. 5.6, together with the asymptotic approximation for large IT. Also shown is the 
numerically computed optimal control. 

We compare the results achieved by piecewise policy timestepping to those achieved by the 
direct control formulation. For clarity, the two discretizations used are 


yn-i-l _ yv 

At 


- min L((T”+^ = 0 

q&Qh ^ 


(5.18) 


for the direct control method and (5.16) for the piecewise constant timestepping method. 

We use policy iteration as in [9| to solve the discrete control problem in (5.18). We use a 
positive coefficient discretization [38] with central differencing used as much as possible. For the 


direct control method, and the piecewise constant policy timestepping method, it is straightforward 
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M = 800 
iV = 50 

J = 5 

M = 1600 
N = 100 

J = 8 

M = 3200 
N = 200 

J = 10 

M = 6400 
N = 400 

J = 15 

M = 12800 
N = 800 

J = 20 

Policy 

(a) 

1.5930 

1.5589 

1.5447 

1.5378 

1.5350 

timestepping 

(b) 


-0.0341 

-0.0141 

-0.0069 

-0.0028 


(c) 



2.41 

2.04 

2.45 

Direct 

(a) 

1.5902 

1.5577 

1.5442 

1.5376 


control 

(b) 


-0.0326 

-0.0135 

-0.0066 



(c) 



2.4167 

2.0333 


Fixed 

(a) 

3.4268 

3.4199 

3.4140 

3.4104 

3.4085 

control 

(b) 


-0.0069 

-0.0059 

-0.0036 

-0.0019 

{q = 1.5) 

(c) 



1.1733 

1.6523 

1.8399 


Table 5.5: The mean-variance test case with parameters as in Table |5.5| with no bankruptcy 
and bounded control. Shown are, for the policy timestepping method, the direct control method, 
and for a fixed constant control: (a) the numerical solution Vk = V{Nk, Mk, Jk',Wo,0); (b) the 
increments Vk — Vk-i; (c) the ratios {Vk — Vk-i)/{Vk-i — Vk- 2 ) for Mk = 800-2*^“^, Nk = 50-2*^“^, 

r 

Jk = 5 ■ v2 (except for the fixed control case, where J = 1). 


to verify that the discretization is monotone, consistent and stable [3811271. The results are shown 
in Table 15.51 

In each step of the policy iteration, the maximum (over parameters) of the discretized differential 
operator at any given mesh-point has to be computed. As the discretization (local upwinding based 
on the coefficients) depends on the control parameter in a discontinuous way, this maximum is found 
by discretizing the control and exhaustive search. This makes the complexity of a single policy 
iteration identical to a single timestep of the constant policy timestepping algorithm. Thus, overall, 
the typically observed 4-6 iterations in every timestep translates into a 4-6 factor of increase in 
the CPU cost of the direct control method compared to the piecewise constant policy timestepping 
technique. Due to this increased cost, we do not show the direct control results for the finest level 
(marked ★). 

The refinements were chosen such that at the coarsest level a single separate refinement of the 
spatial mesh, timestep and control mesh gave comparable (empirical) accuracy improvements. This 
ensures that the data test the convergence order in all three discretization parameters. It is clear 
that the achieved accuracy is almost identical for both methods. 

We also include results for the value achieved with a fixed control, q = 1.5. This is the chosen 
upper bound and the optimal value attained in an interval around Wq = 1, see Fig. 5.6 The 


results are distinctly different from those under the optimal control, which shows that the similar 
performance of policy timestepping and direct control is not a result of the control being constant 
near W = 1. The errors for fixed control are purely due to the time and spatial finite difference 
discretization, and are slightly smaller than those observed in the true optimal control problems. 


6 Conclusions 

This article analyzes the piecewise constant policy timestepping method both from a theoretical 
and an applications perspective. Our main result is that if we use different meshes for each constant 
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policy PDE solve, then convergence to the viscosity solution can be proven even if high order (not 
necessarily monotone) interpolation techniques are used. Essentially, this is because we can view 
the piecewise constant policy timestepping method as the solution to a switching system of PDEs, 
where the coupling between the PDEs occurs only in the zeroth order term. However, this generality 
comes at a price: we must include a finite switching cost in the switching system. Convergence to 
the solution of the original HJB PDE occurs only in the limit as the switching cost tends to zero. 
However, our numerical experiments show that good results are obtained for very small (even zero) 
switching costs. 

The general approach we follow also has superficial similarities with the “semi-Lagrangian meth¬ 
ods” (SLM) of [11] and [T8|. They both make use of the fact that for given coefficients (controls), 
it may be easier to construct monotone schemes together with the underlying mesh, especially in 
more than one dimension. If different controls require different meshes, interpolation of the mesh 
solution is needed in every timestep. In the present method this serves to carry out the optimization 
over solutions with different policies. 

The computational results demonstrate that a smaller error is obtained using the high order 
interpolation, compared to linear interpolation. 

In many practical situations, the local optimization problem at each node is determined by 
discretizing the control and using exhaustive search. In this case, our tests show that piecewise 
constant policy timestepping is more efficient than standard direct control methods, as a similar 
level of accuracy is achieved with less computational effort. This is simply due to the fact that 
piecewise constant policy timestepping is unconditionally stable, and does not require a policy 
iteration to solve nonlinear discretized equations. 

The use of piecewise constant policy timestepping can be useful in situations where generic 
monotone schemes are hard to construct, e.g., in multidimensional settings, whose implementation 
we do not consider here and leave for future work. 

Finally, we note that it is straightforward to implement piecewise constant policy timestepping 
in existing linear PDE solution software. Hence these existing algorithms can be easily converted 
to solve nonlinear HJB equations. 


A Proof of Lemma [T] 

We provide here a proof of Lemma 
Proof. By insertion one gets 


|F(x, (/)(x), 


Fh{x,4>{x) +(,D(j){x),D‘^(j){x))\ 


sup Lq{(j) + C) - sup Lq(j) 
q&Qh q&Q 


< 

sup Lq{(t) + f) -supLg((/)-hO 

+ 

sup Lq{(j) + ^) - sup LqCj) 


q&Qh q&Q 


q£Q q£Q 


by the triangle inequality. From Assumption and the compactness of Q, then the supremum of 
Lqcj) is attained, say at q* , and then 

Lq*(l) - rq*( = Lq*{4> + 0 < SUpLq((/) + 0 < SUp Lq(f> + SUp = Lq*(j) + SUp(-rq)^, 

q£Q qeQ q£Q q&Q 
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hence 


sup Lq((/) + ^) - sup LqCf) 
q£Q q£Q 


< ^sup|rq| 
q&Q 


(A.l) 


Now let be the maximizer of -bq((/> + C). We also have by (2.1) that there is G Qh with 
IqJi — < h. By uniform continuity of the coefficients in q on the compact set Q, there exists a 
function wi so that 

II + ll%* “ II + 1^9? “ '^ 1 *, I + lA; “ I < Wi(x, /i) ^ 0, /i ^ 0 


with the usual vector and matrix norms, and hence 


+ 6 - ‘^i(x,/i)max (l, 1^1, |(/)|, \\D(p\\, \\D'^(j)\\) = max(l, |^|)w 2 (x,/i) 

for a suitably defined 0 ) 2 . Using also that Qh C Q, 


sup Lq{(f> + ^) - max(l, |^|)a; 2 (x,/i) < sup Lq{4> + ^) < sup Lq{(f> + ^), 
q&Q q&Qh q&Q 


so that 


sup Lq{(j) + ^) - sup Lq{4> + ^) 
q&Qh q&Q 


(7)2lx hi ^2 

< max(l, |^|)w 2 (x,/i) <Lj 2 (x,/r)+ ^ ^ + y. (A.2) 


From Assumption!^ and noting that the equation coefficients are uniformly continuous in q, it is 
easily shown that the right hand side of (A.2) is locally Lipshitz in x, independent of h. The result 
then follows from ( A.l|) and (A.2), with an appropriate choice of ooh and o;^. 

□ 
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